Map of Infectious Disease Notifications

Continuing my explortion of Notification of Infectious Disease data, I’ve been working on an interactive map of where disease have been reported in England and Wales. I found this link a useful resource https://journocode.com/2016/01/28/your-first-choropleth-map/.

Aquiring and processing data

I was able to download data from the Notification of Infectious Disease website on the distrubtion of cases by local authority (“Statutory notifiable diseases - number of cases reported in week 51 of 2017”). I’ve limited this experimental analysis to one week of data. I did some refomatting in Excel and saved this as a .csv for speed and convinience, before loading it into R. I then explicitly set NA data to 0 (assuming no notifications means no cases). I also removed the string " UA" from where it appeared in local authority names.

library("reshape2")
library("magrittr")
geo_disease_data<-read.csv("geo_disease_data.csv")
geo_disease_data[is.na(geo_disease_data)] <- 0
geo_disease_data$LocalAuthority<-gsub(" UA", "", geo_disease_data$LocalAuthority)

The next task was to get some map data for local authorities. I played around with a few different maps from the Office for National Statistics Open Geography Portal (I’m new to geo-spatial data). This one seemed to be the best:

library("rgdal")
LocalAuthorities<-readOGR("https://opendata.arcgis.com/datasets/686603e943f948acaa13fb5d2b0f1275_4.geojson")

OGR data source with driver: GeoJSON 
Source: "https://opendata.arcgis.com/datasets/686603e943f948acaa13fb5d2b0f1275_4.geojson", layer: "686603e943f948acaa13fb5d2b0f1275_4"
with 380 features
It has 10 fields

I then merged the two data sets by local authority name.

LA_Diseases<-sp::merge(LocalAuthorities, geo_disease_data, by.x="lad16nm", by.y="LocalAuthority")

The map

The next task was to display the data as a map. For this I used leaflet. This is a wrapper for a Javascrip “App”. I’ve put some comments in the code (after the map) to explain what’s going on at each step.

library("leaflet")

diseases <- names(geo_disease_data)[-1] #Get a list of disease names (the first colmn name here is local authority name)

max_cases <- max(melt(geo_disease_data, id.vars="LocalAuthority")$value) #Get the maximum number of cases for any disease

#Set up a colour scale from 0 to the maximum number of diseases
pal <- colorNumeric(
   palette = "viridis",
   domain = c(0,max_cases)
 )

#Initialise the map with the legend and a control to allow users to select a disease
myleaflet<-leaflet(width = "100%") %>% 
  addLegend( pal = pal, 
             values = 0:max_cases,
             title = "Cases",
             opacity = 1) %>%
  addLayersControl(baseGroups=diseases, 
                   position = "bottomleft", 
                   options = layersControlOptions(collapsed = TRUE))

#Add a layer to the map for each disease
for (active_disease in diseases) 
  {
  myleaflet <- myleaflet %>%
    addPolygons(data=LA_Diseases, 
              fillColor=~pal(LA_Diseases[[active_disease]]),
              fillOpacity = 0.8,
              color = "black",
              weight = 1,
              popup = paste(LA_Diseases$lad16nm, LA_Diseases[[active_disease]],"cases"), #This popup shows the local authority name and number of cases
              group = active_disease
              )
  }    

myleaflet<-myleaflet %>% showGroup("Scarlet.fever") #Pick an interesting disease as default

myleaflet #show the map

Evaluation and conclusions

I have a few problems! Performance on my mobile is not great, but aside from this the main issue is that I don’t have data for all of the local authorities. I thought I would try to find out why so I got a list of local authorities from the disease data and another list from the map.

Unique_LA_Map<-data.frame(LocalAuthority=unique(LocalAuthorities$lad16nm) %>% sort())
Unique_LA_Diseases<-data.frame(LocalAuthority=unique(geo_disease_data$LocalAuthority) %>% sort())

Next, I did some joins to see how many of the names were shared between the data sets and how many were unique to each.

In_Map_Not_Diseases<-data.frame(LocalAuthority=
  Unique_LA_Map$LocalAuthority[!(Unique_LA_Map$LocalAuthority %in% Unique_LA_Diseases$LocalAuthority)])
In_Diseases_Not_Map<-data.frame(LocalAuthority=
  Unique_LA_Diseases$LocalAuthority[!(Unique_LA_Diseases$LocalAuthority %in% Unique_LA_Map$LocalAuthority)])
In_Both<-base::merge(Unique_LA_Map,Unique_LA_Diseases, by="LocalAuthority")

I’ve plotted the results on a Venn diagram. Not the most beautiful, but quite informative.

library(VennDiagram)
grid.newpage()
draw.pairwise.venn(area1 = nrow(Unique_LA_Diseases), 
                   area2 = nrow(Unique_LA_Map), 
                   cross.area = nrow(In_Both), 
                   category = c("n in Diseases Data", "n in Map Data"),
                   fill = c("blue", "red")
                   )

(polygon[GRID.polygon.21], polygon[GRID.polygon.22], polygon[GRID.polygon.23], polygon[GRID.polygon.24], text[GRID.text.25], text[GRID.text.26], lines[GRID.lines.27], text[GRID.text.28], text[GRID.text.29], text[GRID.text.30]) 

276 local authorities are in both data sets, 104 are only in the map data and 12 are unique to the disease data. A large part of the explaination of this is that the map includes scottish local authorities, but the disease data is only for England and Wales. Looking at the data frames, there are some instances where the name is stated differently in each set - despite referring to the same authority e.g. “City of Kingston upon Hull” and “Kingston upon Hull, City of”. And there are some like “East Dorset” that are on the map but don’t seem to have a match at all.

My conclusion is that, whilst making maps of infectious disease occurances is a great way of getting insights out of the data, I would need to access some domain knowledge on mapping local authorities to take this much further.

See my website for more examples. Contact me on bob dot turner dot uk at gmail dot com with comments or questions.

Contains public sector information licensed under the Open Government Licence v3.0.